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ABSTRACT 

We study the final state of the gravitational collapse of uniformly rotating supramassive neutron stars 
by axisymmctric simulations in full general relativity. The rotating stars provided as the initial condition 
are marginally stable against quasiradial gravitational collapse and its equatorial radius rotates with the 
Kepler velocity (i.e., the star is at the mass-shedding limit). To model the neutron stars, we adopt the 
polytropic equations of state for a wide range of the polytropic index as n = 2/3, 4/5, 1, 3/2 and 2. 
We follow the formation and evolution of the black holes, and show that irrespective of the value of 
n (2/3 < n < 2), the final state is a Kerr black hole and the disk mass is very small (< 10~ 3 of the 
initial stellar mass). 

Subject headings: black hole physics - relativity - hydrodynamics - stars: rotation 



1. INTRODUCTION 

Neutron stars are in general rotating. Rotation can 
support neutron stars with higher mass than the maxi- 
mum static limit, producing supramassive stars, as defined 
and numerically computed by Cook et al. (1992, 1994a). 
Supramassive neutron stars may be created (i) when neu- 
tron stars accrete gas from a normal binary companion 
(Cook et al. 1994b), (ii) after the merger of binary neu- 
tron stars (Shibata & Uryu 2000, 2002), and (iii) after 
gravitational collapse of massive stellar core. 

Since viscosity drives any equilibrium star to a uniformly 
■ rotating state, stationary neutron stars are believed to be 
uniformly rotating. The final state after the collapse of 
the marginally stable and uniformly rotating supramas- 
sive neutron stars is the subject of this paper. 

Rotating neutron stars with a density higher than a crit- 
ical value are unstable against gravitational collapse. Such 
' critical density is determined using the turning point the- 
orem (Friedman et al. 1988; Cook et al. 1992). The final 
state of the unstable spherical stars in the adiabatic col- 
lapse is a Schwarzschild black hole. On the other hand, in 
the rotating case, it is not trivial: All the fluid elements 
may not collapse to a Kerr black hole, leaving a fraction 
of the mass around the black hole to form disks. The fi- 
nal state after the collapse of rotating stars is one of the 
fundamental questions in general relativistic astrophysics. 

To clarify the final state of the gravitational collapse of 
rotating neutron stars, numerical simulations in full gen- 
eral relativity are the best approach. Two groups have 
already performed the simulations for relativistic collapse 
of rotating stars (Nakamura 1981; Nakamura et al. 1987; 
Stark & Piran 1985; Piran & Stark 1986). However, they 
have not studied the collapse of marginally stable rotat- 
ing neutron stars which are plausible initial conditions for 
the collapse in nature. Probably this is because numeri- 
cal methods for computation of initial data sets describing 
rapidly rotating neutron stars, as well as numerical tools, 
techniques and sufficient computational resources have be- 
come available only quite recently. Over the last 15 years, 
robust numerical techniques for constructing equilibrium 



models of rotating neutron stars in full general relativity 
have been established (Komatsu et al. 1989; Cook et al. 
1992; Salgado et al. 1994; Stergioulas 1998). More re- 
cently, robust methods for the numerical evolution of the 
coupled equations of Einstein's and hydrodynamic equa- 
tions have been also established (e.g., Shibata 1999b, 2003; 
Font 2000; Font et al. 2002; Siebel et al. 2002, 2003). 

In a previous paper (Shibata et al. 2000), we reported 
the first numerical result for the gravitational collapse, 
which was computed by a three-dimensional numerical im- 
plementation in full general relativity. In that paper, we 
adopted the polytropic equation of state with n = 1 where 
n is the polytropic index, and gave a uniformly rotating 
and marginally stable neutron star at a mass-shedding 
limit (at which the equator of a star rotates with the Ke- 
pler velocity) as the initial condition. The total grid num- 
ber in the simulations was only 153 x 77 x 77 for x — y — z 
(we assumed the equatorial plane symmetry and 7r rota- 
tion symmetry) because of the restricted computational 
resources at that time, and as a result, the equatorial ra- 
dius (polar radius) of the neutron star is covered only by 40 
(23) grid points initially. We found that the collapse leads 
to a black hole (we determined the location of the apparent 
horizon) , and indicated that nonaxisymmetric instabilities 
do not turn on during the collapse. However, we were not 
able to determine the final state of the gravitational col- 
lapse because of the insufficient grid resolution. 

Since nonaxisymmetric instabilities are not likely to be 
relevant during the collapse, the simulation should be car- 
ried out under the assumption of the axial symmetry. 
With this restriction, we could significantly improve the 
grid resolution for a given computational resource. Moti- 
vated by this fact, we have constructed a numerical code 
for axisymmetric numerical simulation in full general rel- 
ativity, which has been already completed (Shibata 2000, 
2003). Because of the restriction to the axial symmetry as 
well as progress in computational resources, we can eas- 
ily increase the grid number 3-5 times as large as that 
in the previous three-dimensional simulation (Shibata et 
al. 2000) even in inexpensive personal computers. As 
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a result, we can search for convergent numerical results 
changing the grid number for a wide range with inexpen- 
sive computational cost. In addition, we adopt a high- 
resolution shock-capturing scheme for evolving the hydro- 
dynamic equations (Shibata 2003), which enables us to 
assess whether shocks play an important role during the 
collapse to a black hole. 

In this paper, we present new numerical results for grav- 
itational collapse computed by the new axisymmetric nu- 
merical implementation. The simulations were carried out 
setting marginally stable equilibrium neutron stars as the 
initial condition. We focus on the collapse of uniformly ro- 
tating supramassive neutron stars at mass-shedding limits 
as before. By exploring rotating stars at mass-shedding 
limits, we can clarify the final state of the collapsed objects 
most efficiently. To investigate the effect of the stiffness of 
equations of state, we adopt polytropic equations of state 
with a wide variety of the polytropic index. The state of 
marginally stable stars which is characterized by the com- 
pactness, angular momentum parameter and density dis- 
tribution depends strongly on the equations of state. This 
implies that the final state after the gravitational collapse 
of rotating stars could depend strongly on the equations 
of state, in contrast to the collapse of nonrotating stars 
in which the Schwarzschild black hole with no disks is the 
unique outcome. To classify the type of gravitational col- 
lapse and its final state, a systematic study for a wide 
variety of equations of state is essential. 

In Sec. 2, we briefly describe our formulation, initial 
data, and spatial gauge conditions. In Sec. 3, we present 
numerical results. In Sec. 4, we provide a summary. 
Throughout this paper, we adopt the units G = c — K = 1 
where G, c and K denote the gravitational constant, speed 
of light and polytropic constant, respectively. We use 
Cartesian coordinates x k = (x, y, z) as the spatial coordi- 
nate with r = yjx 2 + y 2 + z 2 ; t denotes coordinate time. 

2. SUMMARY OF FORMULATION AND INITIAL 
CONDITIONS 

We performed hydrodynamic simulations in general rel- 
ativity for the axisymmetric spacetime using the same for- 
mulation as that used in a previous paper (Shibata 2003), 
to which the reader may refer for details and basic equa- 
tions. 

We assume that neutron stars are composed of the in- 
viscid, ideal fluid. Then, the fundamental variables for the 
hydrodynamic equations are: 

p : rest mass density, 
e : specific internal energy, 
P : pressure, 



can be integrated as 



: four velocity, 



v = 



dx l 
~dt 



(1) 



where subscripts i,j,k,--- denote x,y and z, and p the 
spacetime components. As the fundamental variables to 
be evolved in the numerical simulations, we in addition 
define a density /»*(= pau*e 6 ^) (</> is defined below) and 
weighted four- velocity u,(= (1 + e + P/p)ui) from which 
the total rest mass and angular momentum of the system 



J 



= J d 3 xp*, 
= J d 3 xp*u^ 



(2) 
(3) 



General relativistic hydrodynamic equations are solved us- 
ing the so-called high-resolution shock-capturing scheme 
(Shibata 2003; see Font 2002 for a general review of high- 
resolution shock-capturing schemes) with the cylindrical 
coordinates. 

We neglect effects of viscosity and magnetic fields. The 
dissipation and angular momentum transport timescales 
by these effects are much longer than the dynamical 
timescale unless the magnitude of viscosity and magnetic 
fields is extremely large (e.g., Baumgarte ct al. 2000). 
Thus, neglecting them is appropriate assumption. 

The fundamental variables for the geometry are: 

a : lapse function, 
(3 k : shift vector, 

-fij : metric in 3D spatial hypersurface, 



7 = e 



= det(7y), 



lij = 

Kij : extrinsic curvature. 



(4) 



Aij 



As in the series of our papers, we evolve 7^, 
e~ A ^(Kij —JijK k k ) together with the three auxiliary func- 
tions Fi = S^djjik and the the trace of the extrinsic cur- 
vature K k k with a free evolution code (see Shibata & Uryu 
2002 for the latest version of the formulation) . 

The Einstein equations are solved in the Cartesian coor- 
dinates. To impose the axisymmetric boundary condition, 
the so-called Cartoon method is used (Alcubierre et al. 
2001b): Assuming a reflection symmetry with respect to 
the z = plane, we perform simulations using a fixed uni- 
form grid with the size Nx3xNmx — y — z which covers 
a computational domain as0<a:<L,0<z<L and 
~~ A < y < A. Here, N and L are constants and A = L/N. 
For y = ±A, the axisymmetric boundary conditions are 
imposed. 

The slicing conditions are basically the same as those 
adopted in previous papers (Shibata 1999, 2000, 2003; Shi- 
bata & Uryu 2000, 2002); i.e., we impose an approximate 
maximal slice condition (K k k ~ 0). On the other hand, 
we adopt two spatial gauge conditions for the shift vector. 
One is an approximate minimal distortion (AMD) gauge 
condition [Di(d t ^) — where L>i is the covariant deriva- 
tive with respect to 7^] (Shibata 1999) which has been 
used in our previous works. In contrast with previous pa- 
pers (e.g., Shibata et al. 2000), we used the AMD gauge 
condition without modification. The other is a dynamical 
gauge condition (e.g., Alcubierre et al. 2001a; Lindblom 
& Scheel 2003). In the present work, we impose the dy- 
namical gauge condition with the equation 



^ kl (F l +Atd t F l ), 



(5) 



where At denotes a timestep in numerical computation. 
The second term in the right-hand side of Eq. (5) is in- 
troduced to stabilize numerical computation. With this 
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choice, /3 obeys a hyperbolic-type equation (for a suffi- 
ciently small value of At), because the right-hand side of 
the evolution equation for Fi contains vector Laplacian 
terms as (3 k u + f3 l ik /3 (e.g. Shibata & Nakamura 1995; 
Shibata & Uryu 2002). The outstanding merit of this 
gauge condition is that we can save computational time 
significantly, since we do not have to solve elliptic-type 
equations. In the numerical computations, we adopted 
these two spatial gauge conditions, and found that both 
give the (almost) identical numerical results: As in the case 
of the AMD gauge condition, the dynamical gauge enables 
to carry out a longterm stable simulation irrespective of 
the equations of state. Thus, here, we present numerical 
results with the dynamical gauge condition to demonstrate 
its robustness. 

During numerical simulations, violations of the Hamil- 
tonian constraint and conservation of mass and angular 
momentum are monitored as code checks. Several test 
calculations, including stability and collapse of spherical 
and rotating neutron stars, as well as convergence tests, 
have been described in a previous paper (Shibata 2003). 
Formation of a black hole is determined by finding an ap- 
parent horizon. 

To model supramassive neutron stars, we adopted the 
polytropic equations of state of the form 

P = Kp 1+ ±. (6) 

In this paper, we choose n = 2/3, 4/5, 1, 3/2 and 2 to 
systematically study the effects of stiffness of equations of 
state. During the simulations, we use a T-law equation of 
state as 

P = (r - l)pe, (7) 

where T is the adiabatic constant which is set as 1 + 1/n. 
In the absence of shocks, no heat is generated and the col- 
lapse is adiabatic, preserving the polytropic form of the 
equations of state. This implies that the quantity P/ p T 
measures the efficiency of the shock heating. 

As initial conditions, we gave marginally stable and 
uniformly rotating supramassive neutron stars at mass- 
shedding limits in equilibrium states. To induce gravita- 
tional collapse, we initially reduced the pressure (i.e., K) 
uniformly by 0.5% in all the simulations. Whenever we re- 
duce the pressure, we solve the equations for the Hamilto- 
nian and momentum constraints to enforce them at t = 0. 

Marginally stable supramassive neutron stars of poly- 
tropic equations of state with 2/3 < n < 2 have the com- 
pactness 0.06 < M/R < 0.25 (see Table 1). Typical com- 
pactness of neutron stars is considered to be ~ 0.15-0.2 
(Shapiro & Teukolsky 1983; Glendenning 1996). Thus, the 
present choice of n yields plausible models for marginally 
stable supramassive neutron stars. 

Physical units enter the problem only through the poly- 
tropic constant K, which can be chosen arbitrarily or else 
completely scaled out of the problem. Thus, we only dis- 
play the dimensionless quantities which are defined as 

M* = M«K- n/2 , M = MK' n / 2 , R = RK~ n/2 , 

j = jK~ n , p = P K n , n = nK n , (8) 

where M, R, and denote the ADM mass, equatorial cir- 
cumferential radius and the angular velocity. Hereafter, 
we adopt the units of K = 1 so that we will omit the bar. 



In Table 1, we list the rotating stars at mass-shedding 
limits that we picked up as initial conditions in the present 
simulations. All the quantities are scaled to be nondimen- 
sional using the relation described in Eq. (8). Stabil- 
ity of uniformly rotating polytropes with n = 1, 3/2 and 
2 against gravitational collapse has been already studied 
by Cook et al.(1994). Thus, for these polytropic indices, 
we choose the stars close to the marginally stable point 
based on their results. For n — 2/3 and 4/5, we do not 
know the critical point for the stability. As shown by 
Cook et al. (1994), however, for stiff equations of state 
with n <, 1, the stability of the uniformly rotating stars at 
mass-shedding limits changes near a point where the mass 
is maximum. Thus, we choose the stars of nearly maxi- 
mum mass along the sequence of the uniformly rotating 
star at mass-shedding limits. 

The ratio of the kinetic energy to the gravitational bind- 
ing energy for all the stars that we picked up here is much 
smaller than 0.27 which is a widely believed critical value 
for onset of the dynamical bar-mode instability in a uni- 
formly rotating star (Chandrasekhar 1969). Thus, the 
nonaxisymmetric deformation is unlikely to turn on during 
collapse. This justifies that we assume the axial symmetry. 

3. NUMERICAL RESULTS 
3.1. Prediction 

Before presenting numerical results, we here predict the 
plausible outcome of the gravitational collapse. Such pre- 
diction helps to understand the reason that a result ob- 
tained in a numerical simulation should be the output. 

Because of the axial symmetry of the system (and since 
the fluid is assumed to be inviscid), the mass distribution 
as a function of the specific angular momentum as well 
as the baryon rest-mass and angular momentum are con- 
served throughout the evolution of the system. Using this 
fact, we can predict the final state of gravitational collapse 
from the initial condition. 

We define the mass distribution as a function of the spe- 
cific angular momentum according to (e.g., Stark & Piran 
1987) 

M*(j ) = / d 3 xp*, (9) 

Jj<jo 

where j is a value of the specific angular momentum com- 
puted as xu v {= hu v ) and jo denotes a particular value for 
j. In Fig. 1, we show the mass distribution as a function 
of the specific angular momentum M*(j) as a function of 
j/M. 

To predict the final state of the collapse, we assume that 
(i) a black hole is formed after the collapse, (ii) most of 
the mass elements fall into the black hole, and (iii) the 
value of the nondimensional angular momentum parame- 
ter q = J/M 2 of a formed black hole is nearly equal to the 
value of the system. 

Since the value of q of all the stars that we picked up 
here is smaller than unity and no heat source exist in the 
collapsing star, assumption (i) is quite reasonable, (ii) 
and (iii) are also reasonable because the progenitor of the 
collapse is uniformly rotating, so that the effect of cen- 
trifugal force which prevents a fluid element falling into a 
black hole is important only around the low-density outer 
region of the collapsing stars. 
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According to (i)-(iii), we assume that the mass and Kerr 
parameter of the formed black hole are ss M and w Mq. 
Around the Kerr black hole, there is the innermost stable 
circular orbit (ISCO). All the mass elements of circular 
orbits inside the ISCO have to fall into the black hole. 
This implies that a mass element of the specific angular 
momentum which is smaller than the value at the ISCO, 
iiSCOi has to fall into the black hole. 

Assuming that the mass and Kerr parameter of the 
formed Kerr black holes are M and Mq, we computed 
Jisco for the models listed in Table 1 using the formula 
derived by Bardeen et al. (1972). The numerical results 
are described at the last column of Table 1. In all the 
models, jisco is larger than 2.5M. (Note that it is 2\/3M 
for q = 0.) From Fig. 1, we find that the fraction of the 
mass with jisco > 2.5M is approximately zero (less than 
10~ 3 ) irrespective of n. Therefore, the final state is pre- 
dicted to be a Kerr black hole and the disk mass is very 
small (< 10~ 3 of the initial stellar mass) for any value of 
n between 2/3 and 2. 

3.2. Formation and evolution of black holes 

We performed simulations varying N for a wide range 
as 180-480. This grid number is several times larger 
than that in a previous study (Shibata et al. 2000), and 
enables to check the convergence of the numerical solu- 
tions in detail. For n < 3/2, the equatorial radius of 
marginally stable rotating stars are initially covered by 
N/2 grid points. The polar radius is covered by ~ 0.37V 
in this case. For n = 2, equilibrium stars have a more 
centrally-concentrated density configuration than that for 
stiffcr equations of state. To resolve the central region, we 
arrange the grid spacing with which the equatorial radius 
is covered by 5N/6 grid points initially. In this case, the 
polar radius is covered by « N/2 grid points. 

Numerical computation was in part performed on FA- 
COM VPP5000 in the data processing center of National 
Astronomical Observatory of Japan, but most of the sim- 
ulations were carried out using personal computers of 
Pcntium-4 processors, each of which has 2 Gbytes memory 
and 2.8 GHz clock. Even for N — 480, it takes only about 

5 days to finish a job of <~ 40000 time steps on one of these 
computers. 

As a result of the simulations, we found that irrespective 
of the value of n, the collapse proceeds monotonically to 
be a Kerr black hole. During the collapse, shock heating 
is negligible in the central region. Namely, Pj p T remains 
approximately constant. 

In all the simulations, the apparent horizons were de- 
termined in the late phase of the collapse. As reported 
in a previous paper (Shibata 2003), accuracy of numeri- 
cal results, measured by the violation of the Hamiltonian 
constraint, deteriorates monotonically with time, and the 
computations eventually crashed due to the grid stretch- 
ing around the black holes. However, the grid number 
inside the surface of the apparent horizons in this work is 
large enough to resolve the formation, and evolution of the 
black hole for a duration ~ 20M even without black hole 
excision techniques (W. Unruh 1984, unpublished; Seidel 

6 Suen 1992). The duration is in general longer with bet- 
ter grid resolutions, and with the largest grid number, we 
could determine the final state of the collapse approxi- 
mately. However, to carry out a simulation for more than 



20M after the formation of black holes, excision techniques 
are absolutely necessary. 

We have also checked that the mass distribution as a 
function of the specific angular momentum is conserved 
accurately. In Fig. 1(b), we compare the mass distribu- 
tion at t = and at the formation of apparent horizon for 
n = 2 as an example. The figure shows a good conserva- 
tion of it. 

In Fig. 2, we display the square of the mass of the ap- 
parent horizons Mah as a function of time. Here, Mah is 
defined as 

Mah = ^, (10) 

where A denotes the of the apparent horizon (e.g., 
Cook & York 1990). Figure 2 shows that Mah approaches 
an asymptotic value. It is also found that the numerical 
results are convergent with increase of N. 

Together with the evolution of Mj^ H , in Fig. 2, we plot 
the square of the irreducible mass of the event horizon 
M? r for the Kerr black hole of M and J (dotted horizon- 
tal lines) as 

Mf rr = i (M 2 + VM 4 - J 2 ) . (11) 

Here, as M and J, we adopt the total values of the system 
which are computed from the initial data sets. If the fi- 
nal state of the gravitational collapse is a Kerr black hole 
with negligible disk mass, the mass of the apparent hori- 
zon should approach M; rr . Figure 2 clearly shows that 
Mah asymptotically approaches M; rr . Small deviation of 
the asymptotic value of Mah from Mj rr is likely to be a 
numerical error. Indeed, with improvement of the grid res- 
olution, the final value of Mah appears to converge to Mj rr . 
With the best resolution, the magnitude of the numerical 
error in Mah is less than 1%. This implies that the final 
state of the collapse is the Kerr black hole and the frac- 
tion of the disk mass is very small (within the magnitude 
of numerical error ^ 1%). 

To reconfirm this fact, we display the evolution of the 
fraction of the mass located outside a coordinate radius 
r , M*(r > r )/M», for n — 2 (solid curve), 3/2 (dotted 
curve), 1 (dashed curve), and 2/3 (long-dashed curve) in 
Fig. 3. The result for n = 4/5 is essentially the same, so 
that we omit it. Here, tq is chosen as 0.6M, which is 
approximately equal to the asymptotic value of the coor- 
dinate polar-radius of the apparent horizon in the present 
computations. We found that in our gauge, the coordi- 
nate equatorial-radius is slightly larger than the polar- 
radius, so that the fraction of the mass outside the ap- 
parent horizon is slightly smaller than M*(r > r )/M,,. 
Figure 3 shows that before the collapse, most of the fluid 
elements are located outside r — tq, but during the col- 
lapse, M»(r > r )/M» monotonically decreases and ap- 
proaches zero. This confirms the results in Fig. 2. The 
present results also reconfirm the same conclusion reached 
previously in less-resolved simulations for n = 1 (Shibata 
et al. 2000). 

We note that for n = 2, a small amount of mass elements 
appears to be outside the black hole at the termination 
of the simulations. We infer that they will be eventually 
swallowed into the black hole because the specific angular 
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momentum for the fluid element is not large enough to 
form disks around the black hole. 

Recall that we study the collapse of uniformly rotating 
stars of maximum angular velocity, implying that the effect 
of rotation is taken into account most efficiently There- 
fore, we conclude that the final state after the collapse of 
all the marginally stable and uniformly rotating polytropic 
stars with 2/3<n<2isa Kerr black hole and the disk 
mass is < 10~ 3 M*. 

ft should be noted that for smaller value of n, Mah 
reaches M- lvv more quickly. For n = 2/3, the growth 
timescale of Mah from to <~ Mj rr is <~ QM while for 
n = 2, it is larger than 30M. This reflects that the na- 
ture of the collapse depends strongly on the initial density 
configuration which is determined by the stiffness of the 
equations of state. For stiffer equations of state, the den- 
sity of the initial condition distributes rather uniformly. 
Thus, the collapse proceeds coherently. For softer equa- 
tions of state, on the other hand, the initial condition has 
more centrally-concentrated density distribution with low- 
density outer envelops. Thus, the central region collapses 
to a black hole earlier, and then the outer region falls into 
the black hole spending a longer timescale than that for 
stiffer equations of state. 

To illustrate that the identical numerical results were 
obtained in two different spatial gauge conditions, in Fig. 
4, we display evolution of <f> and p* at r = as well as M AH 
for n = 3/2 with N = 360. The solid and dotted curves 
denote the numerical results in the dynamical and AMD 
gauge conditions, respectively. We see that both results 
are in good agreement. 

Finally, we address the following point: The collapse 
of compact stars to black holes is among the most inter- 
esting processes leading to the production of gravitational 
waves. As pointed out by Stark and Piran (1985), quasi- 
normal modes of a black hole would be excited after the 
formation, and as a result, gravitational waves associated 
with such quasinormal-mode oscillations may be emitted. 
It is an interesting subject to clarify how large the am- 
plitude of gravitational waves is. From this motivation, 
we tried to extract gravitational waves in the simulation, 
but we were not able to do because the amplitude is likely 
much smaller than the typical size of numerical noise of our 
present simulation. The reason that the amplitude is very 
small is the following: (i) the collapse coherently proceeds, 
i.e., almost all the fluid elements collapse to form a black 
hole simultaneously. In such case, the excitation of the 
quasinormal modes of a black hole is likely to be weak, be- 
cause the quasinormal modes are excited by perturbations 
struck after formation of a black hole; (ii) the nondimen- 
sional angular momentum parameter q is not very large 
< 0.7. As Stark and Piran showed that a large amount of 
gravitational waves is emitted for q close to 1. 

4. SUMMARY AND DISCUSSION 

We have reported new numerical results of axisymmet- 
ric simulations for the gravitational collapse of rapidly and 



uniformly rotating supramassive neutron stars to black 
holes in full general relativity. The initial conditions for 
the neutron stars are given using polytropic equations of 
state for a wide range of the polytropic index as n = 2/3, 
4/5, 1, 3/2 and 2. The initial state of the rotating stars 
is marginally stable against the quasiradial gravitational 
collapse and at the mass-shedding limit. The hydrody- 
namic simulations were carried out using a high-resolution 
shock-capturing scheme with the T-law equations of state. 
We have demonstrated that irrespective of the value of 
n (2/3 < n < 2), the collapse monotonically proceeds 
with negligible shock heating, and the final state is a Kerr 
black hole with a small fraction of the disk mass. 

As mentioned in Sec. 3.1, the results obtained in this 
paper can be expected from the initial conditions. In the 
same manner, we can expect the final states of the gravi- 
tational collapse for softer equations of state with n > 2. 
With a large value of n <~ 3, we may model an unstable 
massive stellar core at the final stage of stellar evolution 
and a supermassive star of M <; 10 5 M Q . In Fig. 5, we 
show the mass distribution as a function of the specific an- 
gular momentum of the marginally stable and uniformly 
rotating stars at mass-shedding limits for n = 2.5, 2.9, and 
3. The marginally stable stars for these polytropic indices 
have been already determined by Cook et al. (1994) for 
n = 2.5 and 2.9 and Baumgarte and Shapiro (1999) for 
n = 3. The nondimensional angular momentum param- 
eter q is w 0.39, 0.57, and 0.96 for n = 2.5, 2.9 and 3, 
so that jisco/M for a black hole of mass M and angular 
momentum J is w 3.0, 2.8, and 1.8, respectively. From 
Fig. 5, we can expect that the final state after the col- 
lapse for n = 2.5 is a Kerr black hole and only a small 
fraction of the initial stellar elements (~ 10 -3 M„) forms 
the disks. On the other hand, disks of mass of <; 0.01M* 
and ^ 0.1M* are likely to be formed for n — 2.9 and 3, 
respectively. (If the mass of the black holes is smaller than 
M, jisco is also smaller and, hence, the disk mass could 
be larger. This implies that the numerical fraction men- 
tioned here is the minimum value.) The same conclusion 
for n = 3 has been already drawn by Shibata and Shapiro 
(2002) and Shapiro and Shibata (2002) in a more careful 
analysis. 

The reason why disks are formed for n <; 2.9 is simply 
that the marginally stable stars with polytropic equations 
of state of such large value of n have a large equatorial 
radius with R/M ^> 200, and hence the specific angular 
momentum for a certain fraction of the fluid elements is 
large enough to escape from swallowing into a black hole. 
The present study together with the previous one (Shi- 
bata & Shapiro 2002) shows that nature of the collapse of 
rapidly rotating stars to a black hole depends strongly on 
the equations of state in particular for n ~ 3. 

This work was supported by Japanese Monbuka- 
gakusho Grants (Nos. 13740143, 14047207, 15037204 and 
15740142). 
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(a) (b) 

Fig. 1. — (a) The mass distribution as a function of the specific angular momentum for n = 2/3-2 at t = 0. (b) The same as (a) but at 
t = (solid curves) and at the formation of the apparent horizon (filled circles) for n = 2. This is the result for a simulation with N = 480. 
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Fig. 2. — Evolution of the area of the apparent horizon in units of 16-irM 2 (i.e., the square of the apparent horizon mass) for (a) n = 2/3, 
(b) 4/5, (c) 1, (d) 3/2 and (d) 2. For (a)— (d), the solid, dashed and dottcd-dashed curves are results for N = 360, 240, and 180, and for (e), 
the solid and dashed curves are results for TV = 480, 360, and 240. The dotted horizontal lines denote the area of the event horizon for Kerr 
black holes of a given set of J and M for which we adopt the values of the initial conditions. 
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Fig. 3. — Evolution of the total mass outside a coordinate radius ro for n = 2 (solid curve), 3/2 (dotted curve), 1 (dashed curve), and 
2/3 (long-dashed curve), ro is chosen as ~ 0.6M which is approximately equal to the asymptotic value of the coordinate polar-radius of the 
apparent horizon in the present computations. To plot the curves, we choose the numerical results with N = 360. The time is shown in unit 
of M. 
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Fig. 4. — Evolution of <fr and p* at r = as well as Mj^ H for n = 3/2 with N = 360. The solid and dotted curves denote the results by the 
dynamical and AMD gauge conditions, respectively, and they approximately coincide. 
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